eegmp_analysis/tools/[MEG]PLS/MEGPLS_PIPELINE_v2.02b/MEGPLS_COMPONENTS/View4D_BSR.m

function varargout = View4D_BSR(varargin)
% Viewer for 4D PLS BSRimage .BRIK files.
% - Allows for montages across both slice dimensions and time.
% - Also displays timeseries for selected voxel in BSRimg, GrpAvg or Correlations.
%
% Usage: 
%  View4D_BRIK(anat_img, brik_img, cellarray_of_condition_groupaverages)
%
% Example: 
%  View4D_BRIK('anat+tlrc.BRIK', 'functional+tlrc.BRIK', ...
%    {'GPAVG_L1_0.5-55Hz_0.040s_1.240s+tlrc.BRIK'; 'GPAVG_RT_0.5-55Hz_0.040s_1.240s+tlrc.BRIK'})
%
% Written by: Natasa Kovacevic
% Modified in areas marked by ** (Nov, 2013: M. Cheung)
% Last modified: Jan. 15, 2014


% Copyright (C) 2013-2014, Natasa Kovacevic
%
% This file is a part of the MEG & PLS Pipeline (MEGPLS). For more 
% details, see the documentation included with the software package.
%
% MEGPLS is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License version 2 as published by 
% the Free Software Foundation. This program is distributed in the hope 
% that it will be useful, but WITHOUT ANY WARRANTY; without even the 
% implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License along
% with this program. If not, you can download the license here: 
% <http://www.gnu.org/licenses/old-licenses/gpl-2.0>.


% NK: note use set/get for fields that come prediscribed with objects (e.g. sliders) and for fields that I added to handles. Use setappdata/getappdata for gui
% Begin initialization code - DO NOT EDIT
gui_Singleton = 0;
gui_State = struct('gui_Name',       mfilename, ...
   'gui_Singleton',  gui_Singleton, ...
   'gui_OpeningFcn', @View4D_BSR_OpeningFcn, ...
   'gui_OutputFcn',  @View4D_BSR_OutputFcn, ...
   'gui_LayoutFcn',  [] , ...
   'gui_Callback',   []);
if nargin && ischar(varargin{1})
   gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
   [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
   gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% --- Executes just before View4D_BSR is made visible.
function View4D_BSR_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to View4D_BSR (see VARARGIN)

% Choose default command line output for View4D_BSR
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);



if numel(varargin) < 3 | numel(varargin)>3
   disp('You must supply: anatomical underlay, PLS result bootstrap file and cell array of GPAVG files');
   disp('All files should be in AFNI format');
   disp('Usage: View4D_BSR(anat_img, bsr_img, cellarray_of_condition_groupaverages)');
   disp('Example: View4D_BSR(''anat+tlrc.BRIK'', ''PLS_NR500_1grp_L1vsRT_avg_0.040s_1.240s_LV1+tlrc.BRIK'', {''GPAVG_L1_0.5-55Hz_0.040s_1.240s+tlrc.BRIK''; ''GPAVG_RT_0.5-55Hz_0.040s_1.240s+tlrc.BRIK''})');
   error('wrong number of inputs ');
end

% load anatomical BRIK file
Opt.format = 'vector';
[err,anat,Info,ErrMessage]=BrikLoad(varargin{1},Opt);
handles.anat = anat;

% load PLS BSR image:
Opt.format = 'vector';
[err,img,Info,ErrMessage]=BrikLoad(varargin{2},Opt);
BSR.img = img;
handles.dims = size(BSR.img);  % 4-dimensional (x, y, z, time)
if numel(size(BSR.img))==3, handles.dims(4) = 1; end
handles.BSR = BSR;
handles.Info = Info;
guidata(hObject, handles);

% Sets title to PLS LV info:
%tmp=strfind(Info.HISTORY_NOTE, '***');
% if ~isempty(tmp)
%     tmp=Info.HISTORY_NOTE(tmp:end);
%     loc=strfind(tmp,'\n');
%     finaltmp=tmp(1:(loc(1)-1));
%     set(handles.Title,'String',finaltmp);
% else

    %** Changed to just display full path of BSRimg as title.
    set(handles.Title,'String',varargin{2});
%end

% load condition GPAVG:
for l=1:length(varargin{3})
   Opt.format = 'vector';
   [err,condimg{l},Info,ErrMessage]=BrikLoad(varargin{3}{l},Opt);
end
handles.condimg = condimg;
handles.gpavgname = varargin{3};

%** set initial clamp range & BSR thresholds:
[~, FuncFileName, ~] = fileparts(varargin{2});

if strcmp(FuncFileName(1:7), 'BSRimg_')
    set(handles.TextboxClampRange, 'String', '10');
    set(handles.BSR_posthresh, 'String', '3');
    set(handles.BSR_negthresh, 'String', '-3');
    
elseif strcmp(FuncFileName(1:7), 'SALimg_')
    set(handles.TextboxClampRange, 'String', '0.01');
    set(handles.BSR_posthresh, 'String', '0.001');
    set(handles.BSR_negthresh, 'String', '-0.001');
end


% set view and sliders
x = ceil(handles.dims(1)/2);
y = ceil(handles.dims(2)/2);
z = ceil(handles.dims(3)/2);

set(handles.AxialSlider,'Min',1,'Max',handles.dims(3),'SliderStep',[1/(handles.dims(3)-1) 1/(handles.dims(3)-1)],'Value',z);
set(handles.SagitalSlider,'Min',1,'Max',handles.dims(1),'SliderStep',[1/(handles.dims(1)-1) 1/(handles.dims(1)-1)],'Value',x);
set(handles.CoronalSlider,'Min',1,'Max',handles.dims(2),'SliderStep',[1/(handles.dims(2)-1) 1/(handles.dims(2)-1)],'Value',y);

setappdata(handles.Axialview,'slicenum',z);
setappdata(handles.Sagitalview,'slicenum',x);
setappdata(handles.Coronalview,'slicenum',y);

if numel(size(BSR.img))==3
   handles.dims(4) = 1;
   set(handles.TimeptSlider,'Min',1,'Max',1.001,'SliderStep',[1 1],'Value',1);
else
   set(handles.TimeptSlider,'Min',1,'Max',handles.dims(4),'SliderStep',[1/(handles.dims(4)-1) 1/(handles.dims(4)-1)],'Value',1);
end

% set up montage selection listbox and monatge viewing parameters
set(handles.SelectMontage,'String',{'None','Axial','Sagittal','Coronal'});
set(handles.MontageViewingTable,'data',[1 1 handles.dims(1); 1 1 handles.dims(2); 1 1 handles.dims(3); 1 1 3;]);
% start with first timept
set(handles.Timept,'String','1');

% set up bsr, tbsr, blend
UpdateBlend(handles);

% start using gui
UpdateSlices(x,y,z,handles);
UpdateTimept(1,handles);
guidata(hObject, handles); % Update handles structure

function varargout = View4D_BSR_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;

function UpdateBlend(handles)
% set up bsr, tbsr, blend
% this function is normally called only in the begining and when user
% changes bsr thresholds -> because thhreshold is the only thing that
% alters the blend and the colormap
% clamp bsr image for better viewing

%** get clamp range from new textbox:
bsr_clamp = str2num(get(handles.TextboxClampRange, 'String'));
bsr = handles.BSR.img;
bsr(find(bsr>=bsr_clamp)) = bsr_clamp; bsr(find(bsr<=-bsr_clamp)) = -bsr_clamp;


% get bsr threshold
%** Added check to prevent threshold from exceeding max/min of data or clamp
posthresh = str2num(get(handles.BSR_posthresh,'String'));
negthresh = str2num(get(handles.BSR_negthresh,'String'));

maxvalue = getappdata(handles.BSR_posthresh,'d');
minvalue = getappdata(handles.BSR_posthresh,'c');

if posthresh > maxvalue | posthresh <= minvalue
    posthresh = maxvalue;
    set(handles.BSR_posthresh,'String', posthresh);
end
if negthresh < minvalue | negthresh >= maxvalue
    negthresh = minvalue;
    set(handles.BSR_negthresh,'String', negthresh);
end

% zero out bsr voxels below threshold
tbsr = bsr;
tbsr(find(tbsr<posthresh & tbsr>negthresh)) = 0;
orig_tbsr = tbsr;


% now map tbsr for colormap
d = max(max(max(max(bsr)))); c = min(min(min(min(bsr))));
if c==d, % that would happen only if tbsr is zero everywhere
   tbsr = ones(size(tbsr));
else
   tbsr = round(255*(tbsr-c)/(d-c)+1);
end

[cmap, ignore_pts] = set_colormap(d,c,posthresh, negthresh);
colormap(cmap);

% map grayscale image into midwindow in the colormap
img = handles.anat;
a = min(min(min(img)));
b = max(max(max(img)));
img = round(ignore_pts(1) + (ignore_pts(end)-ignore_pts(1))*(img-a)/(b-a));
img = min(max(img,ignore_pts(1)),ignore_pts(end));
img = repmat(img,[1 1 1 handles.dims(4)]);
blend = zeros([size(tbsr) 3]);

for dim=1:3
   utmp = reshape(cmap(img,dim),size(tbsr));
   otmp = reshape(cmap(tbsr,dim),size(tbsr));
   blend(:,:,:,:,dim)= (orig_tbsr==0).* utmp + (orig_tbsr~=0).* otmp;
end
% attach blend and colormap related infor to be used in other gui
% components
setappdata(handles.BSR_posthresh,'blend',blend);
setappdata(handles.BSR_posthresh,'cmap',cmap);
setappdata(handles.BSR_posthresh,'ignore_pts', ignore_pts);
setappdata(handles.BSR_posthresh,'c',c);
setappdata(handles.BSR_posthresh,'d',d);

function UpdateSlices(x,y,z,handles)
% get current bsr thresholds
posthresh = str2num(get(handles.BSR_posthresh,'String'));
negthresh = str2num(get(handles.BSR_negthresh,'String'));
bsr = handles.BSR.img;
% get  current time point and make sure that it makes sense
tpt = round(get(handles.TimeptSlider,'value'));

% get blend and colormap
blend_allt= getappdata(handles.BSR_posthresh,'blend');
blend = squeeze(blend_allt(:,:,:,tpt,:));
cmap = getappdata(handles.BSR_posthresh,'cmap');
colormap(cmap);

% AFNI IJK coordinates:
set(handles.Coord,'String',[num2str(x) '  ' num2str(y) '  ' num2str(z)]);

% AFNI indexing convention for AFNI viewer & functions starts at 0.
% Note: AFNI reads orientation in voxel-storage order, so LPI in AFNI is RAS (in spatial direction)
% Recall: RAS is orientation of coordinate space of TLRC and MNI templates.
[err, XYZdic] = AFNI_Index2XYZcontinuous([x-1, y-1, z-1], handles.Info, 'LPI');
set(handles.MNI_coord, 'String', [num2str(XYZdic(1)) '  ' num2str(XYZdic(2)) '  ' num2str(XYZdic(3))]);

% udate slider positions
set(handles.AxialSlider,'Value',z);
set(handles.SagitalSlider,'Value',x);
set(handles.CoronalSlider,'Value',y);

setappdata(handles.Axialview,'slicenum',z);
setappdata(handles.Sagitalview,'slicenum',x);
setappdata(handles.Coronalview,'slicenum',y);

% update slice views
col = [0 0.8 0];
axes(handles.Axialview); % make this current axes
tmp=permute(squeeze(blend(:,:,z,:)),[2 1 3]);
image(tmp,'ButtonDownFcn',{@click_axial,handles});
set(handles.Axialview,'ydir','normal','xtick',[],'ytick',[]);
line([x x],[1 handles.dims(2)],'color',[0 1 0],'visible','on');
line([1 handles.dims(1)],[y y],'color',[0 1 0],'visible','on');
text(2,2,['z=' num2str(z)],'fontsize',12,'fontweight','bold','color',[0 1 0.3]);

axes(handles.Sagitalview); % make this current axes
image(permute(squeeze(blend(x,:,:,:)),[2 1 3]),'ButtonDownFcn',{@click_sagital,handles});
set(handles.Sagitalview,'ydir','normal','xtick',[],'ytick',[]);
line([y y],[1 handles.dims(3)],'color',[0 1 0],'visible','on');
line([1 handles.dims(2)],[z z],'color',[0 1 0],'visible','on');
text(2,2,['x=' num2str(x)],'color','g','fontsize',12,'fontweight','bold');
text(0,handles.dims(3)+2,num2str(bsr(x,y,z,tpt)),'fontsize',12,'fontweight','bold','color','c');

axes(handles.Coronalview); % make this current axes
image(permute(squeeze(blend(:,y,:,:)),[2 1 3]),'ButtonDownFcn',{@click_coronal,handles});
set(handles.Coronalview,'ydir','normal','xtick',[],'ytick',[]);
line([x x],[1 handles.dims(3)],'color',[0 1 0],'visible','on');
line([1 handles.dims(1)],[z z],'color',[0 1 0],'visible','on');
text(2,2,['y=' num2str(y)],'color','g','fontsize',12,'fontweight','bold');
text(0,handles.dims(3)+2,num2str(bsr(x,y,z,tpt)),'fontsize',12,'fontweight','bold','color','m');


% update BSR colorbar
axes(handles.colorbar_view);
imagesc([1:256]');  colormap(cmap);
ignore_pts = getappdata(handles.BSR_posthresh,'ignore_pts');

%** Commented out: Modified below to be more robust accomodating different datasets.
%set(gca, 'ydir', 'normal', 'ytick',[1 ignore_pts(1) ignore_pts(end) 256], 'YTickLabel', {num2str(getappdata(handles.BSR_posthresh,'c')), num2str(negthresh), num2str(posthresh), num2str(getappdata(handles.BSR_posthresh,'d'))}, 'XTickLabel', []);
if ignore_pts(1) <= 1
    set(gca, 'ydir', 'normal', 'ytick',[1 ignore_pts(end) 256], 'YTickLabel', {num2str(getappdata(handles.BSR_posthresh,'c')), num2str(posthresh), num2str(getappdata(handles.BSR_posthresh,'d'))}, 'XTickLabel', []);
elseif ignore_pts(end) >= 256
    set(gca, 'ydir', 'normal', 'ytick',[1 ignore_pts(1) 256], 'YTickLabel', {num2str(getappdata(handles.BSR_posthresh,'c')), num2str(negthresh), num2str(getappdata(handles.BSR_posthresh,'d'))}, 'XTickLabel', []);
elseif ignore_pts(1) <= 1 && ignore_pts(end) >= 256
    set(gca, 'ydir', 'normal', 'ytick',[1 256], 'YTickLabel', {num2str(getappdata(handles.BSR_posthresh,'c')), num2str(getappdata(handles.BSR_posthresh,'d'))}, 'XTickLabel', []);
else
    set(gca, 'ydir', 'normal', 'ytick',[1 ignore_pts(1) ignore_pts(end) 256], 'YTickLabel', {num2str(getappdata(handles.BSR_posthresh,'c')), num2str(negthresh), num2str(posthresh), num2str(getappdata(handles.BSR_posthresh,'d'))}, 'XTickLabel', []);
end


% update BSR time series plot
axes(handles.BSR_timeseries);
Tseries = squeeze(handles.BSR.img(x,y,z,:));
if numel(Tseries) > 1
   at = min(Tseries); bt = max(Tseries);
   num_tpts = numel(Tseries);
   tpts = [1:num_tpts];
   plot(tpts,Tseries,'b');
   hold on;
   if bt > at
      e = (bt-at)/20;
      set(gca,'ylim',[at-e bt+e],'xlim',[0 numel(Tseries)+1]);
   else
      set(gca,'ylim',[-1 1],'xlim',[0 numel(Tseries)+1]);
   end
end
% horizontal dotted red lines indicate bsr thresholds
h=line([1 num_tpts], [posthresh posthresh]);
set(h, 'LineStyle', '--', 'Color', 'r');
h=line([1 num_tpts], [negthresh negthresh]);
set(h, 'LineStyle', '--', 'Color', 'r');
hold off;
title(['BSR = ',num2str(Tseries(tpt))],'color','r','fontsize',14);
minval = min(min(Tseries),negthresh); maxval = max(max(Tseries),posthresh);
E = (maxval - minval)/10;
set(gca,'ylim',[minval-E maxval+E]);
% vertical gray line indicates current tpt
h = line([tpt tpt], [minval-E maxval+E]);
set(h, 'LineStyle', '-', 'Color', [0.6 0.6 0.6]);
set(gca,'ButtonDownFcn',{@click_BSR_ts,handles},'xlim',[0 num_tpts+1]);
setappdata(handles.PrintBSRButton,'bsr_ts',Tseries);

% update group-average condition timeseries plots:
axes(handles.GPAVG_timeseries);
num_GPAVG = length(handles.condimg);
tmparray = zeros(num_GPAVG, length(Tseries));
for l=1:num_GPAVG
   tmparray(l,:) = squeeze(handles.condimg{l}(x,y,z,:));
end
tmparray = tmparray';
plot(tmparray); hold on;

for f = 1:length(handles.gpavgname)
    [~, File, ~]  = fileparts(handles.gpavgname{f});
    LegendName{f} = File;
end
clickableLegend(LegendName,'interpreter','none');

YLim = get(gca, 'YLim');
BSRcircleheight = (YLim(2)-((YLim(2)-YLim(1))/10));
circles = find(Tseries > posthresh | Tseries < negthresh);
scatter(circles, BSRcircleheight*ones(size(circles)),'ro');
% vertical gray line indicates current tpt
h = line([tpt tpt], YLim);
set(h, 'LineStyle', '-', 'Color', [0.6 0.6 0.6]);
hold off;
set(gca,'ButtonDownFcn',{@click_GPAVG_ts,handles},'xlim',[0 num_tpts+1]);
setappdata(handles.PrintGPButton,'gpavg_ts',tmparray);

function UpdateTimept(tpt,handles)
% update views
x = getappdata(handles.Sagitalview,'slicenum');
y = getappdata(handles.Coronalview,'slicenum');
z = getappdata(handles.Axialview,'slicenum');
UpdateSlices(x,y,z,handles);

function AxialSlider_Callback(hObject, eventdata, handles)
z = round(get(hObject, 'Value'));
x = round(get(handles.SagitalSlider,'Value'));
y = round(get(handles.CoronalSlider,'Value'));
UpdateSlices(x,y,z,handles);

function SagitalSlider_Callback(hObject, eventdata, handles)
x = round(get(hObject, 'Value'));
z = round(get(handles.AxialSlider,'Value'));
y = round(get(handles.CoronalSlider,'Value'));
UpdateSlices(x,y,z,handles);

function CoronalSlider_Callback(hObject, eventdata, handles)
y = round(get(hObject, 'Value'));
z = round(get(handles.AxialSlider,'Value'));
x = round(get(handles.SagitalSlider,'Value'));
UpdateSlices(x,y,z,handles);

function Timept_Callback(hObject, eventdata, handles)
tpt = round(str2num(get(hObject,'String')));
if tpt<1 || tpt>handles.dims(4)
   prompt={
      ['Invalid Time.'];
      ['Time must be # between: 1 & ',num2str(handles.dims(4))]};
   warndlg(prompt,'Warning:','modal');
end
set(handles.TimeptSlider,'Value',tpt);
UpdateTimept(tpt,handles);

function TimeptSlider_Callback(hObject, eventdata, handles)
tpt = round(get(hObject, 'Value'));
set(handles.Timept,'String',num2str(tpt));
UpdateTimept(tpt,handles);

function click_axial(src,eventdata,handles)%callback for buttonclick
c= ginput(1); 
x = round(c(1,1));
y = round(c(1,2));
z = getappdata(handles.Axialview,'slicenum');
UpdateSlices(x,y,z,handles);

function click_sagital(src,eventdata,handles)
%c = get(gca,'CurrentPoint');
c= ginput(1); 
y = round(c(1,1));
z = round(c(1,2));
x = getappdata(handles.Sagitalview,'slicenum');
UpdateSlices(x,y,z,handles);

function click_coronal(src,eventdata,handles)
%c = get(gca,'CurrentPoint');
c= ginput(1); 
x = round(c(1,1));
z = round(c(1,2));
y = getappdata(handles.Coronalview,'slicenum');
UpdateSlices(x,y,z,handles);

function click_BSR_ts(src,eventdata,handles)
c = ginput(1); 
tpt = round(c(1));
tpt = min(max(1,tpt), handles.dims(4));
set(handles.Timept,'String',num2str(tpt));
set(handles.TimeptSlider,'value',tpt);
UpdateTimept(tpt,handles);

function click_GPAVG_ts(src,eventdata,handles)
c = ginput(1); 
tpt = round(c(1));
tpt = min(max(1,tpt), handles.dims(4));
set(handles.Timept,'String',tpt);
set(handles.TimeptSlider,'value',tpt);
UpdateTimept(tpt,handles);

function Coord_Callback(hObject, eventdata, handles) % user eneters voxel coordinates
coord = round(str2num(get(hObject,'String')));
if (coord(1)<1 || coord(1)>handles.dims(1)) || (coord(2)<1 || coord(2)>handles.dims(2)) || (coord(3)<1 || coord(2)>handles.dims(3))
   prompt={
      ['Invalid Coordinate.'];
      ['X-coord must be # between: 1 & ',num2str(handles.dims(1))];
      ['Y-coord must be # between: 1 & ',num2str(handles.dims(2))];
      ['Z-coord must be # between: 1 & ',num2str(handles.dims(3))]};
   warndlg(prompt,'Warning:','modal');
end
UpdateSlices(coord(1),coord(2),coord(3),handles);

function MNI_coord_Callback(hObject, eventdata, handles)
mni_coord = round(str2num(get(hObject,'String')));
[err, coord] = AFNI_XYZcontinuous2Index(mni_coord, handles.Info, 'LPI', 3);
coord = coord+1; % AFNI indexing convention starts at 0.
UpdateSlices(coord(1),coord(2),coord(3),handles);

function BSR_negthresh_Callback(hObject, eventdata, handles)
thresh = str2num(get(handles.BSR_negthresh,'String'));
if thresh>0, set(handles.BSR_negthresh,'String','0'); end
UpdateBlend(handles);
z = getappdata(handles.Axialview,'slicenum');
x = getappdata(handles.Sagitalview,'slicenum');
y = getappdata(handles.Coronalview,'slicenum');
UpdateSlices(x,y,z,handles);

function BSR_posthresh_Callback(hObject, eventdata, handles)
thresh = str2num(get(handles.BSR_posthresh,'String'));
if thresh<0, set(handles.BSR_posthresh,'String','0'); end
UpdateBlend(handles);
z = getappdata(handles.Axialview,'slicenum');
x = getappdata(handles.Sagitalview,'slicenum');
y = getappdata(handles.Coronalview,'slicenum');
UpdateSlices(x,y,z,handles);

%** Added clamp range textbox
function TextboxClampRange_Callback(hObject, eventdata, handles)
clamp = str2num(get(handles.TextboxClampRange, 'String'));
if clamp<=0, set(handles.TextboxClampRange, 'String', '10'); end
UpdateBlend(handles);
z = getappdata(handles.Axialview,'slicenum');
x = getappdata(handles.Sagitalview,'slicenum');
y = getappdata(handles.Coronalview,'slicenum');
UpdateSlices(x,y,z,handles);

function SelectMontage_Callback(hObject, eventdata, handles)
contents = cellstr(get(hObject,'String')); % get the whole list as a cell array of strings
montage_type = contents{get(hObject,'Value')}; % get current selection
if strmatch(montage_type,{'Axial','Sagittal','Coronal'})
   TableData = get(handles.MontageViewingTable,'Data'); % get User defined parmaeters for the montage
   figure('NumberTitle','off','name',['Montage ' montage_type]); 
   dims = handles.dims;
   % get blend and subsample according to the montage table parameters
   blend = getappdata(handles.BSR_posthresh,'blend');
   subblend = blend(TableData(1,1):TableData(1,2):TableData(1,3),TableData(2,1):TableData(2,2):TableData(2,3),TableData(3,1):TableData(3,2):TableData(3,3),TableData(4,1):TableData(4,2):TableData(4,3),:); % (x y z t rgb)
  
   % determine axes permuation order 
   switch montage_type
      case 'Axial',    ord = [1 2 3 4]; xlab = 'z slice';
      case 'Sagittal', ord = [2 3 1 4]; xlab = 'x slice';
      case 'Coronal',  ord = [1 3 2 4]; xlab = 'y slice';
   end
 
   % permute subblend in awy that is good for the display
   subblend = permute(subblend,[ord(2) ord(4) ord(1) ord(3) 5]); % 5th dim is rgb
   for d=1:4
      span{d} = TableData(ord(d),1):TableData(ord(d),2):TableData(ord(d),3);
   end
      
   % determine vertical image size (dim1) and horizontal image size (dim2)
   dim1 = numel(span{2}) * numel(span{4}); % rows = ysize*time
   dim2 = numel(span{1}) * numel(span{3}); % columns = xsize*zsize
   
   % reshape subblend and show
   subblend = reshape(subblend,dim1,dim2,3);
   image(subblend,'ButtonDownFcn',{@click_Montage,handles});
   set(gca,'ydir','normal','xtick',([1:numel(span{3})]-0.5)*numel(span{1}),'xticklabel',cellstr(int2str(span{3}'))','ytick',([1:numel(span{4})]-0.5)*numel(span{2}),'yticklabel',cellstr(int2str(span{4}'))');
   ylabel('Time'); 
   xlabel(xlab);
   
   % get current gui (x y z t) coordinate
   z = getappdata(handles.Axialview,'slicenum');
   x = getappdata(handles.Sagitalview,'slicenum');
   y = getappdata(handles.Coronalview,'slicenum');
   t = round(get(handles.TimeptSlider,'value'));
   coord = [x y z t];
   
   % tranaslate this coordinate into this image space
   for d=1:4, subcoord(d) = round((coord(d)-TableData(d,1))/TableData(d,2))+1; end
   subcoord = subcoord(ord);
  
   % plot crosshair and 
   horpos = numel(span{2})*(subcoord(4)-1) + subcoord(2);
   verpos = numel(span{1})*(subcoord(3)-1) + subcoord(1);
   horline = line([1 dim2],[horpos horpos],'color',[0 1 0],'visible','on');
   verline = line([verpos verpos],[1 dim1],'color',[0 1 0],'visible','on');
   
   % record current montage window settings
   setappdata(handles.MontageViewingTable,'span',span);
   setappdata(handles.MontageViewingTable,'ord',ord);
   setappdata(handles.MontageViewingTable,'subcoord',subcoord);
   setappdata(handles.MontageViewingTable,'fig',gcf);
   setappdata(handles.MontageViewingTable,'horline',horline);
   setappdata(handles.MontageViewingTable,'verline',verline);
   setappdata(handles.MontageViewingTable,'horpos',horpos);
   setappdata(handles.MontageViewingTable,'verpos',verpos);
end

function click_Montage(src,eventdata,handles)%callback for buttonclick
% % when user click on the montage figure reverse the caluclation from SelectMontage_Callbac
% % to obtain [x y z t] in the original image space
% c = ginput(1);
% horpos = round(c(1,1));
% verpos = round(c(1,2));
% rspan = getappdata(handles.MontageViewingTable,'span');
% ord = getappdata(handles.MontageViewingTable,'ord');
% TableData = get(handles.MontageViewingTable,'Data'); % get User defined parmaeters for the montage
% revord(ord) = [1:4]; % reverse order; this statement defines new variable revord
% for d=1:4, span{d} = rspan{revord(d)}; end
% subcoord(2) = mod(horpos,numel(span{2})); 
% subcoord(4) = (horpos - subcoord(2))/numel(span{2})+1;
% subcoord(1) = mod(verpos,numel(span{1})); 
% subcoord(3) = (verpos - subcoord(1))/numel(span{1})+1;
% revord(ord) = [1:4]; % reverse order; this statement defines new variable revord
% subcoord 
% ord 
% revord
% %subcoord = subcoord(revord)
% for d=1:4, 
%    coord(d) = round( (subcoord(d)-1)*TableData(d,2) + TableData(d,1) ); 
%    coord(d) = min(max(1,coord(d)),handles.dims(d)); 
% end
% set(handles.Timept,'String',num2str(coord(4)));
% set(handles.TimeptSlider,'value',coord(4));
% 
% % redraw crosshairs
% dim1 = numel(span{2}) * numel(span{4}); % rows = ysize*time
% dim2 = numel(span{1}) * numel(span{3}); % columns = xsize*zsize
% horline = getappdata(handles.MontageViewingTable,'horline');
% verline = getappdata(handles.MontageViewingTable,'verline');
% set(horline,'visible','off');set(verline,'visible','off');
% figure(getappdata(handles.MontageViewingTable,'fig'));
% horline = line([1 dim2],[horpos horpos],'color',[0 1 0],'visible','on');
% verline = line([verpos verpos],[1 dim1],'color',[0 1 0],'visible','on');
% setappdata(handles.MontageViewingTable,'horline',horline);
% setappdata(handles.MontageViewingTable,'verline',verline);
% setappdata(handles.Sagitalview,'slicenum',coord(1));
% setappdata(handles.Coronalview,'slicenum',coord(2));
% setappdata(handles.Axialview,'slicenum',coord(3));
% set(handles.Timept,'String',num2str(coord(4)));
% set(handles.TimeptSlider,'value',coord(4));
% UpdateTimept(coord(4),handles);


function PrintGPButton_Callback(hObject, eventdata, handles)
val = get(hObject,'Value');
if val==1
   gpavg_ts = getappdata(handles.PrintGPButton,'gpavg_ts')
   set(hObject,'Value',0);
end

function PrintBSRButton_Callback(hObject, eventdata, handles)
val = get(hObject,'Value');
if val==1
   bsr_ts = getappdata(handles.PrintBSRButton,'bsr_ts')
   set(hObject,'Value',0);
end

function MontageViewingTable_CellEditCallback(hObject, eventdata, handles)
val = round(str2num(eventdata.EditData)); % get cell content as entered by user
inds = eventdata.Indices; % find (i,j) position in the table
% if new value is outside of alowable dims, then change back to previous
% cell content
if ~(val >=1  & val <= handles.dims(inds(1)))
   TableData = get(handles.MontageViewingTable,'Data');
   TableData(inds(1),inds(2)) = eventdata.PreviousData;
   set(handles.MontageViewingTable,'Data',TableData);
end

function MontageViewingTable_CellSelectionCallback(hObject, eventdata, handles)
% hObject    handle to MontageViewingTable (see GCBO)
% eventdata  structure with the following fields (see UITABLE)
%	Indices: row and column indices of the cell(s) currently selecteds
% handles    structure with handles and user data (see GUIDATA)
% NK: note we must have this function even if it doesn't do anything!!!!!





function [cmap,ignore_pts] = set_colormap(max_bsr, min_bsr, pos_bsr_thresh, neg_bsr_thresh)
% written by Jimmy - NK made some changes
%   set the display colormap based on the max/min display values and 
%   the threshold setting.
%
%   The upper colors are coming from the entries of [140:239] of the 
%   255 jet colormap, and the lower colors are from the entries of 
%   [1:100] of the colormap.
%

   range_interval = max_bsr - min_bsr;
   upper_interval = max_bsr - pos_bsr_thresh;
   lower_interval = neg_bsr_thresh - min_bsr; % abs(min_bsr) - abs(neg_bsr_thresh);
   
   %  colormap entries for the upper range values, using the
   %  entries of [140:239] from the 255 jet colormap
   %
   num_upper_colors = 0;
   if (upper_interval > 0)
      num_upper_colors = round(upper_interval / range_interval * 255); 
      cmap_size = round(255 * num_upper_colors/100);
      first_color_idx = round(140 / 255 * cmap_size);
      last_color_idx = first_color_idx + num_upper_colors - 1;
      uppermap = jet(cmap_size);
      upper_colors = uppermap(first_color_idx:last_color_idx,:);
   end;
   
   %  colormap entries for the lower range values, using the 
   %  entries of [1:100] from the 255 jet colormap
   %
   num_lower_colors = 0;
   if (lower_interval > 0)  
      num_lower_colors = round(lower_interval / range_interval * 255); 
      cmap_size = round(255 * num_lower_colors/100);
      first_color_idx = 1;
      last_color_idx = num_lower_colors;
      lowermap = jet(cmap_size);
      lower_colors = lowermap(first_color_idx:last_color_idx,:);
   end;
   
   cmap = zeros(256,3);
   cmap(1:255,:) = jet(255);
  
   ignore_pts = [(num_lower_colors+1):(255-num_upper_colors)];
   
   if (num_lower_colors > 0),
     cmap(1:num_lower_colors,:) = lower_colors;
   end;
   if (num_upper_colors > 0),
     if ~isempty(ignore_pts)
       cmap((ignore_pts(end)+1):255,:) = upper_colors;
     else
       cmap((255-num_upper_colors+1):255,:) = upper_colors;
     end
   end
   if ~isempty(ignore_pts)
     %cmap(ignore_pts,:) = ones(length(ignore_pts),3) * 140/255;
     cmap(ignore_pts,:) = gray(numel(ignore_pts)); 
   end
   cmap(256,:) = [1 1 1];






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% Create functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function BSR_posthresh_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
   set(hObject,'BackgroundColor','white');
end

function BSR_negthresh_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
   set(hObject,'BackgroundColor','white');
end

function MNI_coord_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
   set(hObject,'BackgroundColor','white');
end

function Coord_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
   set(hObject,'BackgroundColor','white');
end

function TimeptSlider_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
   set(hObject,'BackgroundColor',[.9 .9 .9]);
end

function Timept_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
   set(hObject,'BackgroundColor','white');
end

function CoronalSlider_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
   set(hObject,'BackgroundColor',[.9 .9 .9]);
end

function SagitalSlider_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
   set(hObject,'BackgroundColor',[.9 .9 .9]);
end

function AxialSlider_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
   set(hObject,'BackgroundColor',[.9 .9 .9]);
end

function SelectMontage_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


%** Added clamp range textbox
% --- Executes during object creation, after setting all properties.
function TextboxClampRange_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end